/*
http://www.matthiasbook.de/papers/parallelfractals/mandelbrot.html
http://www.math.utah.edu/~alfeld/math/mandelbrot/mandelbrot.html#algorithm
http://glinda.lrsm.upenn.edu/~weeks/software/mand.html
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mytime.h"

#define MAXITER 100000
#define THRESOLD 2.0
#define npoints 4092

float uni (void);
void rinit (int);

struct complex{
    double creal;
    double cimag;
};

int min(int a, int b) {
    return((a < b) ? a : b); 
}


int main(int argc, char **argv) {
    int i, j, numinside;
    double area_seq, error_seq;
    double area_llc, error_llc;
    double area_mpi, error_mpi;
    struct complex z, c[npoints];

    double ztemp;
    int numoutside, gnumoutside;
    int nt;                             /* No. of threads */
	int MPI_NUMPROCESSORS, MPI_NAME;
    CLOCK_TYPE chrono;
	double t_llc, t_seq, t_mpi;

/*
 *  1. Generate npoints random points in the complex plane
 */

    LLC_printMaster ("\n\n\n*************** NUMPROCESSORS = %d ***************************\n\n", LLC_NUMPROCESSORS);
    rinit (54321);
    for (i=0; i<npoints; i++) {
        c[i].creal = -2.0+2.5*uni();
        c[i].cimag = 1.125*uni();
    }

/*
 *  2. Monte Carlo sampling 
 *
 *    2a. Outer loop runs over npoints, initilaise z=c
 *
 *    2b. Inner loop has the iteration z=z*z+c, and threshold test
 */
	CLOCK_Start(chrono);
	
	numoutside = 0;
    for(i = 0; i<npoints; i++) {
        z.creal = c[i].creal;
        z.cimag = c[i].cimag;
        for (j = 0; j < MAXITER; j++) {
            ztemp = (z.creal * z.creal) - (z.cimag * z.cimag) + c[i].creal;
            z.cimag = z.creal * z.cimag * 2 + c[i].cimag;
            z.creal = ztemp;
            if (z.creal * z.creal + z.cimag * z.cimag > THRESOLD) {
                numoutside++;
                break;
            }
        }
    }
    numinside = npoints - numoutside;

/*
 *  3. Calculate area and error and output the Results
 */

    area_seq = 2.0 * 2.5 * 1.125 * numinside / npoints;
    error_seq = area_seq / sqrt(npoints);
	CLOCK_End(chrono, t_seq);
/*
 *  4. Parallel (with llc)
 */

	CLOCK_Start(chrono);
 
	numoutside = 0;
    #pragma omp target device(mpi)
    #pragma omp parallel for default(none) reduction(+:numoutside) private(i,j,ztemp,z) shared(nt,c)
    for(i = 0; i<npoints; i++) {
        z.creal = c[i].creal;
        z.cimag = c[i].cimag;
        for (j = 0; j < MAXITER; j++) {
            ztemp = (z.creal * z.creal) - (z.cimag * z.cimag) + c[i].creal;
            z.cimag = z.creal * z.cimag * 2 + c[i].cimag;
            z.creal = ztemp;
            if (z.creal * z.creal + z.cimag * z.cimag > THRESOLD) {
                numoutside++;
                break;
            } 
        } /* for j */
    } /* for i */
  
	numinside = npoints - numoutside;

/* *  5. PARALLEL llc: Calculate area and error */
    area_llc = 2.0 * 2.5 * 1.125 * numinside / npoints;
    error_llc = area_llc / sqrt(npoints);
	
	CLOCK_End(chrono, t_llc);

/*
 *  6. Parallel (with MPI)
 */

    MPI_Comm_size(MPI_COMM_WORLD, &MPI_NUMPROCESSORS);
	MPI_Comm_rank(MPI_COMM_WORLD, &MPI_NAME);
	
	CLOCK_Start(chrono);
 
	numoutside = 0;
    for(i = MPI_NAME; i<npoints; i+= MPI_NUMPROCESSORS) {
        z.creal = c[i].creal;
        z.cimag = c[i].cimag;
        for (j = 0; j < MAXITER; j++) {
            ztemp = (z.creal * z.creal) - (z.cimag * z.cimag) + c[i].creal;
            z.cimag = z.creal * z.cimag * 2 + c[i].cimag;
            z.creal = ztemp;
            if (z.creal * z.creal + z.cimag * z.cimag > THRESOLD) {
                numoutside++;
                break;
            } 
        } /* for j */
    } /* for i */

    MPI_Allreduce(&numoutside, &gnumoutside, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);        
	numinside = npoints - gnumoutside;

/* *  7. PARALLEL MPI: Calculate area and error */
    area_mpi = 2.0 * 2.5 * 1.125 * numinside / npoints;
    error_mpi = area_mpi / sqrt(npoints);
	
	CLOCK_End(chrono, t_mpi);

/* *  8. RESULTS */

    printf("(%d/%d)SEQ: Area of Mandlebrot (npoints = %d) set = %16.12f +/- %16.12f Time %g sec.\n", 
            LLC_NAME, LLC_NUMPROCESSORS, npoints, area_seq, error_seq, t_seq);
    printf("(%d/%d)LLC: Area of Mandlebrot (npoints = %d) set = %16.12f +/- %16.12f Time %g sec.\n", 
            LLC_NAME, LLC_NUMPROCESSORS, npoints, area_llc, error_llc, t_llc);
    printf("(%d/%d)MPI: Area of Mandlebrot (npoints = %d) set = %16.12f +/- %16.12f Time %g sec.\n",
            LLC_NAME, LLC_NUMPROCESSORS, npoints, area_mpi, error_mpi, t_mpi);

    LLC_printMaster ("%d\t%g\t#llc_plot0 MANDEL: N = %ld. [seq_time(%g)/llc_time(%g)]\n", 
            LLC_NUMPROCESSORS, (t_seq/t_llc), npoints, t_seq, t_llc);
    LLC_printMaster ("%d\t%g\t#llc_plot1 MANDEL: N = %ld. [seq_time(%g)/mpi_time(%g)]\n", 
            LLC_NUMPROCESSORS, (t_seq/t_mpi), npoints, t_seq, t_mpi);
    LLC_printMaster ("%d\t%g\t#llc_plot2 MANDEL: N = %ld. [mpi_time(%g)/llc_time(%g)]\n", 
            LLC_NUMPROCESSORS, (t_mpi/t_llc), npoints, t_mpi, t_llc);

}
